home *** CD-ROM | disk | FTP | other *** search
/ Mac Magazin/MacEasy 21 / Mac Magazin and MacEasy Magazine CD - Issue 21.iso / Wissenschaft & Technik / yorick12vr1-nofpu folder / include / demo1.i < prev    next >
Text File  |  1995-07-26  |  8KB  |  294 lines
  1. /*
  2.    DEMO1.I
  3.    1-D hydro code written in Yorick language
  4.  
  5.    $Id: demo1.i,v 1.1 1993/08/27 18:50:06 munro Exp $
  6.  */
  7. /*    Copyright (c) 1994.  The Regents of the University of California.
  8.                     All rights reserved.  */
  9.  
  10. /* ------------- Set a few parameters --------------------------- */
  11.  
  12. /* Set initial density and temperature.  */
  13. rho0= 1.1845e-3;             /* g/cc dry air at NTP */
  14. RT0= 298.16;                       /* Kelvin at NTP */
  15.  
  16. /* Set number of zones for hydro calculation, and
  17.    total length of column.  */
  18. n_zones= 200;
  19. column_length= 100.0/*cm*/;
  20.  
  21. /* A conversion factor.  */
  22. R_gas= 8.314e7;    /* erg/K/mole ideal gas constant */
  23.  
  24. /* Gamma-law gas equation of state parameters.  */
  25. gamma= 1.4;                 /* Cp/Cv for air is 7/5 */
  26. gammaM1= gamma - 1;
  27. A_bar= 1.1845e-3/*g/cc dry air*/ * 24.5e3/*cc/mole*/;
  28. K_to_v2= 1.e-6*R_gas/A_bar;     /* (cm/ms)^2/Kelvin
  29.                                    for dry air      */
  30. /* Local temperature units are (cm/ms)^2. --
  31.    Note dependence of temperature units on A_bar.  */
  32. RT0*= K_to_v2;
  33.  
  34. /* ------------- Set initial conditions --------------------------- */
  35.  
  36. /* Set initial conditions for hydro calculation.  */
  37. func reset
  38. {
  39.   extern RT, z, v, M, p;
  40.  
  41.   RT= array(RT0, n_zones);           /* temperature */
  42.   z= span(0, column_length, n_zones+1);     /* zone
  43.                                boundary coordinates */
  44.   v= array(0.0, dimsof(z)); /* zone bndy velocities */
  45.  
  46.   /* The column consists initially of n_zones zones
  47.      of equal size containing equal amounts of gas. */
  48.   M= rho0*abs(z(2)-z(1));         /* mass/area/zone */
  49.  
  50.   /* The pressure array includes a pressure before and
  51.      after the column, p(0) and p(n_zones+1), which
  52.      will be used to set pressure boundary conditions.
  53.      (Velocity boundary conditions will be applied by
  54.      setting v(0) and v(n_zones).)
  55.      Note that the units of this pressure are
  56.           g*(cm/ms)^2/cc.  */
  57.   p= array(rho0*RT0, n_zones+2);
  58. }
  59.  
  60. /* ------------- Set boundary conditions --------------------------- */
  61.  
  62. func sound
  63. /* DOCUMENT sound
  64.      Set up the initial conditions for evolve to launch a weak sound wave.
  65.  
  66.    SEE ALSO: shock, evolve
  67.  */
  68. {
  69.   extern bc0_v, bc0_time, bc0_p, bc0_Z;
  70.   bc0_v= 0.1*sin(span(0, 2*pi, 100));
  71.   bc0_time= span(0, 1.0, 100);
  72.   bc0_p= bc0_Z= [];
  73.   reset;
  74. }
  75.  
  76. func shock
  77. /* DOCUMENT sound
  78.      Set up the initial conditions for evolve to launch a strong wave, which
  79.      steepens into a shock as it propagates.
  80.  
  81.    SEE ALSO: sound, evolve
  82.  */
  83. {
  84.   extern bc0_v, bc0_time, bc0_p, bc0_Z;
  85.   bc0_v= 10.0*sin(span(0, 2*pi, 100));
  86.   bc0_time= span(0, 1.0, 100);
  87.   bc0_p= bc0_Z= [];
  88.   reset;
  89. }
  90.  
  91. sound;
  92.  
  93. func nobc {
  94.   extern bcN_v, bcN_p, bcN_Z;
  95.   bcN_v= bcN_p= bcN_Z= [];
  96. }
  97.  
  98. func hardbc {
  99.   extern bcN_v, bcN_p, bcN_Z;
  100.   bcN_v= 0;
  101.   bcN_p= bcN_Z= [];
  102. }
  103.  
  104. func softbc {
  105.   extern bcN_v, bcN_p, bcN_Z;
  106.   bcN_p= rho0*RT0;
  107.   bcN_v= bcN_Z= [];
  108. }
  109.  
  110. func matchbc {
  111.   extern bcN_v, bcN_p, bcN_Z;
  112.   softbc;
  113.   bcN_Z= rho0*sqrt((gammaM1+1)*RT0);
  114. }
  115.  
  116. /* ------------- Define the main function --------------------------- */
  117. /* The DOCUMENT comment will be printed in response to:  help, evolve */
  118.  
  119. func evolve(time1, time0)
  120. /* DOCUMENT evolve, time1
  121.          or evolve, time1, time0
  122.      Step the hydro calculation forward to TIME1,
  123.      starting with the initial conditions in the
  124.      RT, z, and v arrays at time TIME0 (default 0.0
  125.      if omitted).  The calculation also depends on
  126.      the constants M (mass/area/zone) and gammaM1
  127.      (gamma-1 for the gamma-law equation of state).
  128.      The pressure array p is updated in addition to
  129.      the primary state arrays RT, z, and v.
  130.  
  131.      Boundary conditions are specified by setting
  132.      either a boundary pressure or a boundary
  133.      velocity at each end of the fluid column.
  134.      bc0_v   - Boundary velocity at z(0), or []
  135.                if z(0) has pressure BC.
  136.      bc0_p   - Boundary pressure beyond z(0).
  137.      bc0_time  - If bc0_v or bc0_p is an array,
  138.                  bc0_time is an array of the same
  139.                  length specifying the corresponding
  140.                  times for time dependent BCs.
  141.      bc0_Z   - Acoustic impedance at z(0) if bc0_v
  142.                is nil (default is 0).
  143.      bcN_v, bcN_p, bcN_time, and bcN_Z have the same
  144.      meanings for the z(n_zones) boundary.
  145.  
  146.      The worker routines OutputResults and
  147.      TakeStep must be supplied.
  148.  */
  149. {
  150.   if (is_void(time0)) time0= 0.0;
  151.  
  152.   for (time=time0 ; ; time+=dt) {
  153.     dt= GetTimeStep();
  154.     SetBoundaryValues, time, dt;
  155.     OutputResults, time, dt;
  156.     if (time >= time1) break;
  157.     TakeStep, dt;
  158.   }
  159. }
  160.  
  161. /* ----------------- compute time step --------------------- */
  162.  
  163. func GetTimeStep(dummy)
  164. {
  165.   dz= abs(z(dif));
  166.   dv= abs(v(dif));
  167.   cs= sqrt((gammaM1+1)*RT);
  168.   return min( dz / max(courant*cs, accuracy*dv) ) * dt_multiplier;
  169. }
  170. /* Set reasonable default values for courant and
  171.    accuracy parameters.  */
  172. courant= 2.0;  /* number of cycles for sound signal
  173.                   to cross one zone -- must be >=2.0
  174.                   for numerical stability */
  175. accuracy= 3.0; /* number of cycles for zone volume to
  176.                   change by a factor of ~2 -- must be
  177.                   >1.0 to avoid possible collapse to
  178.                   zero volume */
  179. dt_multiplier= 1.0;
  180.  
  181. /* ----------------- set boundary conditions ------------------ */
  182.  
  183. func SetBoundaryValues(time, dt)
  184. {
  185.   vtime= time + 0.5*dt;  /* velocity is 1/2 step
  186.                             ahead of pressure, z */
  187.  
  188.   /* boundary at z(0) */
  189.   if (!is_void(bc0_v)) {
  190.     /* velocity BC */
  191.     v(1)= BCinterp(bc0_v, bc0_time, vtime);
  192.  
  193.   } else if (!is_void(bc0_p)) {
  194.     /* pressure BC */
  195.     p(1)= BCinterp(bc0_p, bc0_time, time);
  196.     /* acoustic impedance BC */
  197.     if (!is_void(bc0_Z))
  198.       p(1)-= sign(z(2)-z(1))*bc0_Z*v(1);
  199.   }
  200.  
  201.   /* boundary at z(n_zones) (written here as z(0)) */
  202.   if (!is_void(bcN_v)) {
  203.     v(0)= BCinterp(bcN_v, bcN_time, vtime);
  204.  
  205.   } else if (!is_void(bcN_p)) {
  206.     p(0)= BCinterp(bcN_p, bcN_time, time);
  207.     if (!is_void(bcN_Z))
  208.       p(0)+= sign(z(0)-z(-1))*bcN_Z*v(0);
  209.   }
  210. }
  211.  
  212. func BCinterp(values, times, time)
  213. {
  214.   if (numberof(times)<2) return values(1);
  215.   else return interp(values, times, time);
  216. }
  217.  
  218. /* ----------------- produce output ----------------------- */
  219.  
  220. func OutputVPlot(time, dt)
  221. {
  222.   extern cycle_number;
  223.   if (time==time0) cycle_number= 0;
  224.   else cycle_number++;
  225.   if (!(cycle_number%output_period)) {
  226.     fma;
  227.     plg, v, z;
  228.     zx= sqrt((gammaM1+1)*RT0)*time;
  229.     if (zx>max(z)) {
  230.       /* make a dotted "ruler" at the location a sound wave
  231.      would have reached */
  232.       zx= 2*max(z)-zx;
  233.       if (zx<min(z)) zx= 2*min(zx)-zx;
  234.     }
  235.     pldj, zx, min(bc0_v), zx, max(bc0_v), type="dot";
  236.   }
  237. }
  238. output_period= 8;  /* about 50 frames for wave to
  239.                       transit 200 zones */
  240.  
  241. /* OutputResults can be switched among several possibilities */
  242. OutputResults= OutputVPlot;
  243.  
  244. /* ----------------- 1-D hydro worker ----------------------- */
  245.  
  246. func TakeStep1(dt)
  247. {
  248.   /* velocities 1/2 step ahead of coordinates */
  249.   z+= dt * v;
  250.   dv= v(dif);
  251.   dz= z(dif);
  252.  
  253.   /* Compute artificial viscosity.  */
  254.   q= array(0.0, n_zones+2);
  255.   q(2:-1)= q_multiplier * M * max(-dv/dz, 0.0)*abs(dv);
  256.  
  257.   /* Apply 1st law of thermodynamics to compute
  258.      temperature change from work p*v(dif) done
  259.      on zone.  Note that p is not time-centered
  260.      properly here.  */
  261.   RT-= (dt * gammaM1/M) * (p+q)(2:-1) * dv;
  262.  
  263.   /* Update the pressures from the updated densities
  264.      M/z(dif) and temperatures RT.  Note that p(0)
  265.      and p(-1) updated by SetBoundaryValues.  */
  266.   p(2:-1)= M * RT/dz;
  267.  
  268.   /* Apply Newton's 2nd law to update velocities.  */
  269.   v-= (dt/M) * (p+q)(dif);
  270. }
  271. q_multiplier= 1.0;
  272.  
  273. TakeStep= TakeStep1;
  274.  
  275. /* ----------------- simple interface ----------------------- */
  276.  
  277. func demo1
  278. /* DOCUMENT demo1
  279.      run the 1-D hydrocode demonstration.  An X window should pop up
  280.      in which a movie of a wave propagating down a shock tube appears.
  281.      Use the 'sound' and 'shock' commands to set two different interesting
  282.      initial conditions; the default is 'sound'.
  283.  
  284.    SEE ALSO: sound, shock, evolve
  285.  */
  286. {
  287.   window, 0, wait=1;
  288.   limits;
  289.   reset;
  290.   evolve, 5;
  291. }
  292.  
  293. /* ----------------- end of demo1.i ----------------------- */
  294.